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Abstract. We consider a two-dimensional lattice model of equilibrium statistical mechanics, 
using nearest neighbor interactions based on the matching conditions for an aperiodic set of 16 
Wang tiles. This model has uncountably many ground state configurations, all of which are 
nonperiodic. The question addressed in this paper is whether nonperiodicity persists at low but 
positive temperature. We present arguments, mostly numerical, that this is indeed the case. In 
particular, we define an appropriate order parameter, prove that it is identically zero at high 
temperatures, and show by Monte Carlo simulation that it is nonzero at low temperatures. 



1. Introduction 

Certain aUoys are beheved to exhibit, at low temperature, a state of thermal equilibrium 
which is solid but not crystalline (as determined for instance by X-ray diffraction), a state 
called quasicrystalline [1]. 

In part to explain the observed diffraction patterns, it is common to model the energy 
ground state of such a material by using (aperiodic) tilings [1, 2, 3]. But ever since such 
models have been proposed [4, 5] there has been the need to determine whether they 
actually are useful to explain the behavior of materials at positive temperature, that is, it 
is unclear whether such models exhibit a phase transition from the usual disorder at high 
temperature to an aperiodically ordered state at low but positive temperature. 

We will analyze a two-dimensional lattice "tiling" model with appropriate energy 
ground states in an attempt to make some progress in this problem. The lattice model is 
of a standard form [6, 7]. In fact the specific model has been discussed already, by Leuzzi 
and Parisi [8], though they concentrated on the degeneracy of the energy ground state 
of the model, and a possible connection with the nonequilibrium behavior of glasses and 
similar materials. See also [9]. 

We will not be giving a proof of a phase transition in this model. Indeed, there are very 
few models for which one can prove phase transitions. A proof of a quasicrystalline phase 
may well require a new basic technique. But at least it would be useful to have a good 
order parameter, and to develop some intuition about the nature of the order, on which 
to base a future argument for a transition. This is our goal. We will introduce such an 
order parameter for our model, prove that it vanishes identically at high temperature, and 
present numerical evidence that it is nonzero below some critical value of the temperature. 
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2. The model and main results 



2.1. An overview 

We consider a model of interacting particles on the lattice Z^, where each site can be in 
one of 16 different states. The possible one-particle states are identified with Ammann's 
16 "prototiles" depicted in Fig. 1 below. Each of these prototiles is a unit square in M^, 
centered at the origin, whose edges each carry a label from one of the sets S = {1,2} or 
L = {3, 4, 5, 6}. (Rotations are not allowed.) Notice that the horizontal edges of a prototile 
are either both of type S, or both of type L. Similarly for the vertical edges. 



1 

2 1 
2 



4 

1 1 

6 



4 

1 2 

5 



6 

1 2 
3 



5 

2 2 
3 



1 

6 4 
1 



2 

5 4 
1 



2 

3 6 
1 



2 

3 5 
2 



3 

5 3 
5 



3 

5 5 
6 



5 

6 3 
5 



5 

6 5 
6 



5 

6 6 
4 



6 

4 5 
6 



6 

4 6 
4 



Fig. 1. Ammann's 16 prototiles. 



Denoting Ammann's set of prototiles by A, a configuration of our particle system is given 
by a map o" : j e i— > cr^ e A. Such a configuration a can also be represented by a 
collection of tiles {j + crj : j e 7?}, where j + aj is a translated copy of the prototile 
CTj . This collection defines a covering of by squares (with labeled edges) whose centers 
lie on Z^, and whose interiors are mutually disjoint. In such a covering, any pair of tiles 
{i + ai,j + aj} that share a common edge assigns two labels to this edge. If the two labels 
disagree, then the set {i,j} is referred to as a defect of a. A configuration or covering 
without defects will be called an ^-tiling, or tiling for short. 

It is shown in [10] that such tilings exist, and that they are all nonperiodic. The 
nonperiodicity follows from the existence of an infiation rule, which is related to the sub- 
stitution rule S ^ L and L — > 5'L for Fibonacci sequences. (Fibonacci sequences are 
2-sided- infinite words, obtained from the sequence 

SL^ LSL SLLSL LSLSLLSL SLLSLLSLSLLSL ^ ... 

via translations and limits.) More precisely, all tiles in a given row of an ^-tiling have 
the same vertical type [S or L), and their horizontal type defines a Fibonacci sequence. 
Similarly for the columns. Furthermore, the inflation map is invertible on the set Q of 
all A-tilings. Besides nonperiodicity, this also implies that Q carries a unique translation- 
invariant probability measure A [3]. 

Since all ^-tilings are nonperiodic, no finite "patch" determines a tiling uniquely. 
However, large patches of a tiling are strongly correlated, even if they are arbitrarily far 
apart. The question that motivated our analysis is whether this "long range order" is still 
present in a thermodynamic ensemble, where typical tile configurations have a small but 
positive density of defects. A thermodynamic ensemble is a probability measure on the 
space S of tile configurations a : Z'^ ^ A, indexed by an inverse temperature /? > 0. A 
more precise definition will be given below. Roughly speaking, our measure v/j assigns a 
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relative weight exp[—l3H{a)] to a tile configuration a, where H{a) is the number of defects 
of a. The minimizers of H, also referred to as ground state configurations, are precisely 
the A-tilings. (It is not difficult to prove that any state obtained as a limit (3 ^ oo gives 
full measure to these configurations [11, 12].) 

To test for long range order, we only consider the type {S or L) of the four tile edges. 
This identification by type defines an equivalence relation "~" on A, with four equivalence 
classes. Given a configuration cr, an A-tiling 7, and a finite region A C Z^, denote by 
(/)A(cr, 7) the fraction of sites j G A where CTj ~ 7j. The limit as A | Z^, if it exists, will 
be denoted by (^(a, 7), and we define V'(cr) = 0(a, 7)(iA(7). One of our goals is to study 
the quantity 

QM7)= / %^rfz.^(a), 7e^, (2.1) 

Js 

which measures how much a typical tile configuration at inverse temperature (3 aligns with 
the ground state configuration 7. In the absence of any preference, Qp is the constant 
function 1, or equivalently, the "order parameter" 

QiP)= fQph)HQph)]dXh) (2.2) 
Jg 

is identically zero. We will prove that this is the case for sufficiently small (3 > 0. This 
result is independent of the choice of boundary condition in defining through the energy 
function H. 




Our remaining results are purely numerical. In what follows, we choose as boundary 
conditions (at infinity) a fixed A-tiling r. It is important to keep this in mind. Fig. 2 
shows the values of q{f3) obtained via Monte Carlo simulations, for tile configurations of 
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size N X N, with A^, a power of 2, ranging from 32 to 256. These results suggest that q{P) 
becomes negative as (3 is increased past a critical value (3c ~ 2.4. If correct, this would 
imply the existence of a phase transition, from a disordered state for /3 < /3c to an ordered 
state for P > Pc , where translation invariance is broken. As one would expect, Qp{'j) 
takes its maximum at 7 = r. Fig. 3 shows the simulated values of Qp{r). 

Similar signs of an order-disorder transition were found in [8], using Monte Carlo 
simulations for 8 < < 32, with free (but eventually frozen) boundary conditions. At 
these values of A^, there is evidence that the phase transition is of second order, with a 
power-law or logarithmic divergence of the specific heat (as /3 — > /3c), depending on the 
model used to fit the data. Our numerical results for 32 < < 256 clearly favor the 
second alternative, if either. In fact, we find a slowdown in the increase of the specific 
heat, suggesting that the phase transition is of third or higher order. Our simulated values 
for the energy per tile, and for the specific heat, are shown in Figs. 4 and 5. The curve 
labeled dE/dT is the (discrete) derivative of the energy, for A^ = 256, and the other three 
curves in Fig. 5 were obtained from the energy fiuctuations. 
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Fig. 4 The average energy per tile. Fig. 5 The specific heat. 



Fig. 6 shows a two-point correlation Cp{d) as a function of the (vertical) separation 
d. For /3 > /3c , the correlation approach a nonzero constant that depends on /3, indicat- 
ing again the existence of long range order at low temperature. The correlation length 
C^"'^(l/4) is shown in Fig. 7. It diverges roughly like (/3c — /3)~^ as /3 t /3c ■ Such power 
law behavior is again similar to what is observed in models with ordered low temperature 
phases, except that the exponent 7 is unusually large. By contrast, the correlations for 
/3 ~ /3c decay roughly like exp(— cd^/^), at least in the observed range. This suggests 
that there is no renormalization fixed point (nontrivial scaling limit) associated with the 
critical value /3c . Instead, there seems to be a "line of fixed points" for /3 > /3c . This is a 
feature that is better known in models with a continuous (internal) symmetry, such as the 
Xy-model. 
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Fig. 6 Correlation: Cp{d) versus d. Fig. 7 Correlation length ""^(1/4). 



2.2. Further observations 

We now give a description of the ground states, that can be used to compute and visuahze 
the function Qp . It should also provide some insight into the behavior of the model at low 
temperature. 

Without the normalizing factor -i/'"^ in the integral (2.1), our function Qp is analogous 
to the order parameters used in models that have periodic ground states and/or a compact 
internal symmetry group. By being a function on Q, it is implicitly covariant under any 
symmetry of the model, including translations. If none of the symmetries are broken, 
then the corresponding "entropy" q{l3) vanishes. Thus, Qp seems to be a natural order 
parameter in models with a large number of ground states. For the tiling model considered 
here, this number is uncountable. Nevertheless, Qp is easy to compute. The reason is 
that, as we will see below, can be identified with the torus T^, where translations 

{T^ a)i = ai^j act by irrational rotations. Each equivalence class [a] = {x G A : x ~ a} of 
A corresponds to one of four disjoint rectangles in a covering of T^, and 'ilj~^{a)(f)f^{a, 7) is 
the probability that the rectangle [aj] contains the point T^'j, for a randomly chosen site 
j e A. 

To be more precise, given any X,Y E {S, L}, denote hy X x Y the set (equivalence 
class) of all prototiles in A whose horizontal and vertical edges are of type X and Y, 
respectively. Let R = {SxS, SxL, LxS, LxL}. Then to every tile configuration a : — > A 
we can associate a function [a] : I? ^ R, by setting [a]j = [aj]. If 7 is an A-tiling, then 
[7] is a product of two Fibonacci sequences /c t— > and /c 1-^ , in the sense that 
[tIj — ^ji ^ yj2 ■ Conversely, any product of two Fibonacci sequences can be obtained in 
this way. Thus, the set of all such i?-tilings [7] is the product T x T , where T denotes 
the set of all Fibonacci sequences. Since (5/3(7) only depends on the equivalence class [7], 
it suffices to find a convenient description for the Fibonacci sequences. (In fact, [7] = {7} 
for almost every tiling 7, but we will not use this here.) 
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One such description is the foUowing [2, p. 128]. Given a real number -d, and a partition 
{J{S), J{L)} of the circle T = M/Z, we can associate with any angle a e T a sequence 
X e {S, L}^ by setting 

_fS, ifa + MeJ{S); . . 

'"'~\L, ila + MeJiL). ^^-"^^ 

Let now be the inverse golden mean, 1? = 1(^-1). The Fibonacci sequences are 
obtained by choosing either J{L) = [0, and ^(5") = ["J?,!), or else J{L) = (0,1?] and 
J{S) — {'&, 1]. The two sets of sequences differ only by a countable set, corresponding to 
angles a + k-d that are zero (modulo 1). This set has measure zero, so we can ignore it. 
Thus, for our purposes, can be identified with the circle T. 

In this representation, translations {T^x)m = Xm+k on JF become irrational rotations 
R^a = a + k'd on the circle. Similarly, Z^-translations on x are represented by 
irrational rotations R^a — {ai + oli + j2^) on the torus T^. In both cases, the unique 
invariant measure is Lebesgue measure. We will describe later how these properties can 
be used for numerical computations. 

While a ground state configuration 7 determines a point a. on the torus, an n x n 
"patch" of 7 determines a rectangular neighborhood of a of linear size 0{n~^\ Thus, a 
typical low temperature configuration determines locally a pair of "fuzzy" angles. Unlike 
in the XY-model, the energy associated with a gradual change of the column (row) angle 
over a horizontal (vertical) distance d does not decrease with d. This follows from the fact 
that the density of letter-mismatches between two Fibonacci sequences is asymptotically 
proportional to the angle difference. A similar argument may apply in the other directions, 
based on the "slanted" Fibonacci sequences described in the remark below. Thus, it seems 
plausible that the model can maintain long range order at low temperatures, despite the 
fact that there are uncountably many ground states. 

By analogy with the Xy- model, one might ask about the existence of vortices and/or 
dipoles in our tiling model. The question is meaningful only at reasonably low temperature, 
since the angles are ill defined at high temperature. So isolated vortices are unlikely to 
play a major role. But in a model with slowly varying angles, a dipole can be associated 
with two successive crossings through a fixed value. Such dipoles (horizontal and vertical) 
are a prominent feature at temperatures near /?c . But they arc rarely isolated and thus 
hard to analyze systematically. Besides these horizontal/vertical dipoles, one can also 
observe "slanted" dipoles whose ends are single defects (not type mismatches), and whose 
connecting line has a slope near . These slanted dipoles seem to be the main source of 
entropy at very low temperatures. Their density has no visible singularity over the range 
of temperatures considered, so they do not appear to play a major role in the observed 
phase transition. 

Remark. In a different representation of the Ammann tilings [10], the 16 prototiles are 
rectangles, whose L-edges and S'-edges have lengths "i? and 1 — t?, respectively. In this rep- 
resentation, it is possible to replace the numeric edge-labels by two types of line segments, 
say blue and green, transverse to the edges, such that a perfect tiling is characterized by 
the blue (and similarly the green) segments combining into a parallel sequence of straight 
lines, known as Ammann bars. The bars are slanted, with slopes ■Q^^ ^ and it should not 
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be too surprising that the spacings between the blue (as well as the green) bars define 
a Fibonacci sequence. The same can be done with square tiles, except that the bars are 
only piecewise linear. This shows that the matching rules for the Ammann tiles enforce, 
primarily and in a direct way, products of Fibonacci sequences. 

2.3. The model 

In this section, we give a more detailed description of the model and show that the order 
parameter (2.2) vanishes for small f3 > 0. We recall that our simulations were carried out 
with boundary conditions given by a tiling r. For the purpose of this section, r could be 
any configuration in S. Thus, we shall suppress the dependence on r in our notation. 

We start by considering finite regions k <zl? . The energy RpXp^ of a configuration 
cr G 5 is defined as the number of defects of a that intersect A. Given a real number /? > 0, 
and a finite subset A of Z^, the Gibbs state (for A) at temperature with boundary 
condition r, is the functional that assigns to a continuous function / (for the product 
topology) on S the value 

(/)/3,A = E • (2.4) 

CTe<SA 

Here, 5a is the set of configurations cr e 5 that agree with r outside A, and .^/3,a is a 
normalization constant, determined by the condition (1)/3,a = 1- Taking a limit K\ 1? 
along squares defines a Gibbs measure on 5, 

= |im(/)^,A. (2.5) 

By well known results in the theory of lattice models [13, 14], the measure for small 
positive /? is translation invariant, ergodic, and mixing. In particular, does not depend 
on the choice of boundary condition r. (For large the limit may have to be taken along 
subsequences, and it can depend on r.) 

Let V be any translation-invariant probability measure on 5, and consider the space 
f2 = 5 X ^, equipped with the product measure ijl = This measure is invariant under 

translations T^{a,^) = {T^a,T^^). Thus by the ergodic theorem [15], if is any function 
in L^(f2), then the orbit averages 

^ n— 1 n — 1 

<^n(^'7)=4^ E E <^(^'(^'7)) (2.6) 

ji=-n j2 = -n 

converge //-almost everywhere to a function in L^{Q), as n — > oo. In what follows, let 
(^(a, 7) = 9{(Tq ~ 7o), where 6'(true) = 1 and ^(false) = 0. Now assume that u is mixing. 
Then, by a standard result in ergodic theory [16, p. 228], the measure fi is ergodic. As a 
result, 4> is the constant function with value (j)dfx. 

This shows that the function defined in (2.1) is well defined, as long as the Gibbs 
measure i/jj is translation invariant. Furthermore, if vp is mixing, then this function is 
identically 1, and q{P) = 0. As mentioned above, this holds for sufficiently small /? > 0. 
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We note that these arguments do not show that is well defined for all values of 

although this seems likely to be true. (It is true for instance at zero temperature, as we 
will see later.) In any case, if translation invariance is broken at some /3 > 0, then a phase 
transition has to occur. 

The correlations described earlier are given by 

r/3(z) = {VoVi)f3 - {Vo)f3{Vi)p , ieZ\ (2.7) 

where Vi = y o T*, and where F : <S — > R is defined as follows. Denote by Xj{a) the 
horizontal type (L or S) of the tile aj in a configuration a. Define f{S) = —1 and 
/(L) = 1. Then V{a) is the value of f{xj{a)), averaged over all sites j in some fixed finite 
region containing the origin. The data in Fig. 6 are only for the correlation Cp{d) = 
in the vertical direction i — {0,d). The correlations for V in other directions have been 
computed as well, but they are not shown here. They have large oscillations (related to 
the Fibonacci sequence) and are less convenient for estimating a correlation length. We 
also considered a purely probabilistic measure of correlations, namely the relative entropy 
of the random variable (j, a) i-^ Xj{a) and its translates. The results are not qualitatively 
different from those found via (2.7). 



3. Computations 

Using the correspondence between ground states (modulo equivalence) and points on the 
torus, the overlap (2.6) of a configuration a with a tiling 7 can now be written as follows. 
Let a be the point on defined by 7. Writing [aj] — Xj x yj , we have 

^ n— 1 n—1 

4>n{cr,l) = ^ xiJ{xj),ai+ji'd)x{J{yj),a2+j2i^), (3.1) 

3i = -n j2 = -n 

where b 1— > x{B, h) denotes the indicator function of a set S C T. The integral of 0„ over 
q; e is given by 

^ n—1 n—1 

V'n(^) = ^ E E \J{^MJ{y3)\- (3.2) 

Ji = -n i2 = -n 

As was shown earlier, ip^^cj)^, — > 1 as — > 00, for small (3. The limit can also be 
computed at zero temperature. In this case, the double sum in (3.1) factorizes into a 
product of simple sums. Without loss of generality (due to translation invariance), we 
can assume that a corresponds to the torus point 0. By using the ergodicity of irrational 
rotations, one finds that (/)^ — > a.e. on ^ x ^, with ^(cr, 7) = (p{ai)(p{a2)i where is the 
piecewise linear function 

m^|l-2|t|, if|t|<l-^?; 
^ [ 2i? — 1 , otherwise. 

The functions V'n converge a.e. to the constant k^, where k = + (1 — d)'^. 
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We expect that Qpi^y) ~ ip{ai)(fi{a2) for large values of /3. This is indeed observed 
numerically, but this is to be expected in a finite system. For comparison, we show in 
Fig. 8 the computed values of Qpi^y), as a function of the point a G T^, for the inverse 
temperatures 2.3 and 2.4. (The observed transition for = 256 is between these values.) 




Fig. 8 Overlap Qp with the different ground states. 

We recall that our simulations were carried out with boundary conditions given by a tiling 
r. This is why (and where) the graphs in Fig. 8 have a single maximum. If we had used 
other boundary conditions, then the limit of cr^ as /9 oo would be a mixture of pure 
tiling states, possibly very complicated, and the system might resemble a spin glass, as 
was observed for instance in [8]. In this context, we should mention that the number of 
defect- free configurations on an x lattice square is bounded by e'^^, as was already 
described in [8] . The number of such configurations that can be extended to a full A-tiling 
is only 0{N^). 

In our numerical computations, we evaluate the sum in (3.1) on a finite 256 x 256 grid 
of points a G T^. (Choosing a finer grid gives no significant improvement.) Averaging 
(f)n over a then yields i/'^ . This is done for each individual tile configuration a in a 
collection T,2n,/3, obtained via Monte Carlo simulation. Then we perform the integral (2.1) 
by averaging over the configurations in E2n,/3 • 

Each of our ensembles Sat,/? contains 10'^ configurations. For increased flexibility, 
they were computed beforehand and stored for analysis later [17]. The configurations 
in Sat,/? are separated from each other by at least M^r^^ Monte-Carlo steps (updates of 
individual tiles) , where Mjv,/3 was determined by monitoring overlaps with the appropriate 
starting configuration, and settling times for various observables, to eliminate any visible 
dependence or bias. To give an example, M256,2.45 ~ 2.3 * 10^^. The initial configurations 
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for N = 256 were obtained by slowly cooling a random configuration. Patches of the 
resulting configurations were also used to generate the starting points for N < 256. 
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